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FOREWORD 


This is a final report on the research project, "Analysis and 
Computation of Internal Flow Field in a Scramjet Engine," for the 
period January 1, 1987 to December 31, 1987. Special attention 
during this period was directed to "Investigation of Supersonic 
Chemically Reacting and Radiating Channel Flow." The work was 
supported by the NASA Langley Research Center (Computational 
Methods Branch of the High-Speed Aerodynamics Division) through 
the grant NAG-1-423. The grant was monitored by Drs. Ajay Kumar 
and J. Philip Drummond of HSAD-Computational Methods Branch, the 
work, in part, was supported also by the Old Dominion University's 
ICAM Project through NASA grant NAG-1-363; this grant was 
monitored by Dr. Samuel E. Massenberg, University Affairs Officer, 
NASA Langley Research Center, Hampton, Virginia 23665. 



INVESTIGATION OF SUPERSONIC CHEMICALLY 
REACTING AND RADIATING CHANNEL FLOW 


By 

Mortaza Mani 1 and Surendra N. Tiwari 2 

SUMMARY 

The two-dimensional time dependent Navier-Stokes equations are used to 
investigate supersonic flows undergoing finite rate chemical reaction and 
radiation interaction for a hydrogen-air system. The explicit multi-stage 
finite volume technique of Jameson is used to advance the governing equa- 
tions in time until convergence is achieved. The chemistry source term in 
the species equation is treated implicitly to alleviate the stiffness asso- 
ciated with fast reactions. The multi-dimensional radiative transfer equa- 
tions for a nongray model are provided for a general configuration and then 
reduced for a planer geometry. Both pseudo-gray and nongray models are used 
to represent the absorption-emission characteristics of the participating 
species. 

The supersonic inviscid and viscous, nonreacting flows are solved by 
employing the finite volume technique of Jameson and the unsplit finite 
difference scheme of MacCormack. The specified problem considered is of the 
flow in a channel with a ten degree compression-expansion ramp. The calcu- 
lated results are compared with the results of an upwind scheme. The pro- 
blem of chemically reacting and radiating flows are solved for the flow of 
premixed hydrogen-air through a channel with parallel boundaries, and a 
channel with a compression corner. Results obtained for specific conditions 
indicate that the radiative interaction can have a significant influence on 
the entire flow field. 

^Graduate Research Assistant, Department of Mechanical Engineering and Me- 
chanics, Old Dominion University, Norfolk, Virginia 23529. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23529. 
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Chapter 1 


INTRODUCTION 

In the last several years, there has been a great deal of research 
toward development of a hypersonic transatmospheric vehicle. At the 
NASA Langley Research Center, the hydrogen-fueled supersonic combustion 
ramjet (scramjet) engine has been a strong candidate for propelling such 
a vehicle. Both experimental and numerical techniques are being 
employed for a better understanding of the complex flow field in 
different regions of the engine. Numerical modeling of the flow in 
various sections has proven to be a valuable tool for gaining more 
insight into the complex nature of these flows [1-5] . 

During the past two decades, a tremendous progress has been made in 
the field of radiative energy transfer in nonhomogeneous nongray gaseous 
systems. In recent years, radiation heat transfer has received 
attention because of its application in fire and combustion research, 
entry and reentry phenomena, and hypersonic propulsion. In the 
hypersonic propulsion system, the temperature ranges from 1000-5000 K. 
In this range, various nonsymmetric molecules such as H 2 O, CO 2 and OH 
become highly radiative participating. Infrared absorption and emission 
of thermal radiation is a consequence of coupled vibrational and 
rotational energy transitions. Diatomic molecule is the simplest 


*The numbers in brackets indicate references. 


1 



2 


molecule that will undergo such transitions. However, symmetric 
diatomic molecules, such as H 2 , 0 2 and N 2 , have no permanent dipole 
moment and thus are transparent to infrared radiation. For unsymmetric 
diatomic and tri atomic molecules, such as OH, CO, C0 2 and H 2 0, the 
infrared spectrum will consist of fundamental vibration-rotation bands 
occurring at the fundamental vibrational frequencies of the molecule 
followed by the overtone and combination bands [6]. 

In the past, radiative transfer analysis, due to the complexity of 
the formulations and the computer resource requirements, was limited to 
one-dimensional formulations. Even for one-dimensional cases, the non- 
gray radiative heat transfer calculations required enormous amount of 
computational time. Important works in nongray one-dimensional 
formulation are reviewed in Refs. 6-9. 

Since the late 1960's, efforts have been directed toward 
formulating efficient and accurate multi-dimensional equations for 
radiative transfer. Latko and Pomraning [10] suggested a synthesis 
method for solving two-dimensional time dependent radiative transfer 
equations. The synthesis method is an attempt to reduce the computer 

time requirements by constructing the two-dimensional problem from a 
small number of one-dimensional calculations. The method is alternated 
from cycle to cycle. On odd cycles, the two-dimensional calculation is 
constructed from N one-dimensional x-direction (y=const.) cuts. On even 
cycles, the two-dimensional calculation is constructed from M one- 
dimensional y-direction (x-const.) cuts. Berg and Crosbie [11] 
developed an exact formulation for the radiative flux and emissive power 
for a two-dimensional finite planer, absorbing-emitting gray medium in 
radiative equilibrium. Exact expressions were obtained for a medium 
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subjected to the following types of boundary conditions: (A) cosine 

varying collimated radiation, (B) a strip of collimated radiation, (C) 
cosine varying diffuse radiation, and (D) a uniform temperature strip. 
The solution for the cosine varying collimated radiation model was used 
to construct the solutions for the other boundary conditions. The two- 
dimensional equations were reduced to one-dimensional equations by the 
method of separation of variables. 

Tsai and Chan [12] presented a general formulation of radiative 
heat flux and its divergence for multi -dimensional radiative problems 
involving nongray absorbing-emitting gases. The expressions obtained 

are in terms of total band absorptance rather than the spectral absorp- 
tion coefficients. Thus the expressions are more compact. Modest [13] 
developed a new multi-dimensional model to calculate the spectrally- 
integrated total radiative-flux for a molecular gas band based on the 
solution of two simple differential equations. This model employs the 
exponential-wide band model and, therefore, considerably reduces the 
numerical efforts required. Yuen and Wong [14] solved two-dimensional 
radiative transfer equations for gray medium by a point alocation method 
in which the temperature profile is expressed as a polynomial of a 
successively higher order. It was shown that the technique provides a 
rapid convergence in comparison to the Hottel's Zonal Method [15] with 
the same number of unknowns. Their approach represents a reduction in 
computational time by about one order of magnitude. This technique 
requires the evaluation of finite number of single integrals for a 
complete solution to the problem. 

In the combustion temperature range, some diatomic and triatomic 
molecules are highly radiative participating species. 


Various 
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investigators have studied the effect of radiative transfer for channel 
flows. Martin and Hwuang [16] solved the energy equation for steam 

flowing between two parallel black walls. The flow was assumed to be 
steady and the radiation transfer in the flow direction was neglected. 
It was shown that the radiative flux peaks at a small distance from the 
wall, instead of at the wall. This effect was also noted by Viskanta in 
a gray analysis of radiation and convection between plates of constant 
temperatures [17]. This is because at the wall the effect of positive 
radiant heat flux from the wall is partially cancelled by the negative 
flux from the layers of hot gas next to the wall. At small distances 
into the stream, however, the flux from both the wall and the hot gas 
combine to give a maximum heat flux. Kobiyama et al. [18] studied the 
problem of combined radiation and convection for compressible laminar 
flow between two isothermal parallel plates. A comparison between 
temperature profiles calculated with the treatment of one- and two- 
dimensional radiation shows a considerable temperature difference at the 
entrance region of the heating zone. The problem of combined convection 
and radiation in a rectangular duct was also studied by Im and Ahluwalia 
[19] for compressible turbulent flows. The moment method was employed 
in this study to solve the radiative flux equation. This method reduces 
the general radiation transport equations to a set of equations in x, y, 
and z directions. It was concluded that the radiative transfer causes 
the thermal boundary layer to grow and skin friction to decrease. The 
velocity profile was not affected by the radiative heat transfer. 

Chung and Kim [20] solved a two-dimensional combined mode heat 
transfer problem by using the finite element technique. The effect of 
scattering was also included in the radiative formulation. It was 
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concluded that the standard Galerkin finite elements may be used if the 
convection domination is relatively small (Re Pr < 1000). If conduction 
energy transfer dominates over radiation, there are few effects of 
optical thickness on the temperature profile in the absence of 
scattering. For converging channels, the radiation effect on the 
temperature profiles is insignificant even when conduction and 
convection are small. 

Tiwari and Singh [21, 22] investigated the transient radiative 
interactions of nongray absorbing-emitting species in laminar fully 
developed flows between two parallel plates. The particular species 

considered were OH, CO, C0 2 , H 2 0 and different mixture of these 

species. Their results demonstrated that, H 2 0 is a highly radiation 

participating species as compared to C0 2 , CO and OH. The effects of 
radiation increase with increasing plate spacing, and the radiative 
transfer is more pronounced at higher wall temperatures and pressures. 
It was also shown that optically thin limits overestimate the influence 
of radiation. Soufiani and Taine [23] studied the H 2 0-air mixtures for 
the above geometry and reached the same conclusions. 

James and Edwards [24] added the nongray radiation described by the 
exponential model for molecular gas bands to the numerical solution of 
turbulent combustion of methane in a planner, enclosed, jet-diffusion 

flame. The planner jet of methane was injected with velocity Uf ue i into 
a stream of air flowing with velocity U a1 - r parallel to the fuel. 

Diffusion-controlled combustion occurred in the mixing region of the 
jet. Plane-parallel isothermal black walls symmetrically located above 
and below the jet formed the combustion chamber. A soot-free flame was 
assumed to exist so that the molecular gas bands determined the thermal 
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radiative transfer to the walls. In this study, 40 percent of the 
computation time was devoted to the radiation calculation. The goal of 
the study was to show the effect of radiation upon, combustion 

temperatures for different channel sizes. Three channel widths of 0.2, 
2.0 and 20.0 meters were selected. It was noted that as the channel 
width reduced, the effects of radiation were also reduced, and for the 
channel width of 0.2m, there was no effect on the temperature profile. 

For the numerical investigation of chemically reacting and radiat- 
ing flows, an appropriate chemistry model must also be selected. 

Depending on the ratio of the chemical and fluid dynamic time scales, 
the suitable chemistry model could be a frozen flow model, a finite rate 
model, an equilibrium model, or a complete reaction model. In general, 
the finite-rate model is the most accurate one. In the last several 
years, a number of finite rate chemistry models for hydrogen-air systems 
have been introduced in the literature. Rogers and Schexnayder [25] 
proposed as many as 60 reaction paths in their model; this is certainly 
one of the most complete representatives of hydrogen-air reaction. 
Unfortunately, the enormous number of reaction paths and chemical 
species involved in the model makes it unfeasible for numerical 
investigation of engineering problems. Intermediate level models are 
reduced to 12 species, 25 reaction paths, eight species and eight 

reaction paths [26]. Except for some inaccuracies during the ignition 
delay period, the eight reaction models perform as well as the 25- 
reaction path model. Although these models are less tedious than the 
60-path model, they are expected to be too costly for use in routine 

parametric studies. The global two-step chemistry model of Rogers and 
Chinitz [27] is an inexpensive and attractive model for primary 
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investigation of reacting flows. This node! was deduced by fitting the 
temperature history of a 28-reaction model [25] used in a series of 
constant-pressure stream-tube calculations. There are a number of 
limitations to this model, such as ignition phase inaccuracy and a 

tendency to overpredict the flame temperature. But as pointed out 
earlier, it is considered to be an appropriate model for the initial 
parametric study of overall mixing and extent of combustion. 

The global two-step chemistry model was used successfully by 
several investigators to solve chemically reacting supersonic flows. 
Drummond et al . [2] used the global chemistry model to solve the flow in 
a rapid expansion nozzle. The governing equations describing the flow 
were solved by the two-stage Runge-Kutta method for integrating in time, 
and a Chebyshev spectral method for integrating the equations in 
space. The results were compared with the two finite difference schemes 
of Adam-Moulton and MacCornack. The comparison showed that the spectral 
method with the Runge-Kutta Scheme gives the same accuracy on much 

coarser grids as compared to the finite-difference procedure. This 

results in a significant gain in the computational efficiency. Bussing 
and Murman [28] solved the time-dependent Navier-Stokes equations for 
supersonic reacting flows. Several efficient acceleration techniques 
were used for calculating steady state chemically reacting supersonic 
flows. The techniques included preconditioning the conservation 

equations, and a preconditioned multiple-grid accelerator. Chitsonboon 
et al. [4] used two-dimensional parabolized Navier-Stokes and 

parabolized species equations to investigate supersonic chemically 
reacting flows related to scramjet-engine configurations. A linearized, 
fully-coupled, fully-impl icit finite difference algorithm was used to 
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develop a computer code to solve the governing equations by marching in 
space rather than time. Results obtained by using the parabolized 
formulation were compared with the results obtained by using the 
elliptic equations. The comparisons indicated fairly good agreement 
between the results of the two formulations. 

A more realistic chemistry model was used by Drummond [5] in 
numerical simulation of a supersonic chemically reacting mixing layer. 
To explore the behavior of such flows, detailed physical models of 
convective and diffusive mixing and finite rate chemical reaction in 
supersonic flow were developed. The finite rate chemistry model 
consisted of eighteen reaction paths and nine species. In this study, 
two numerical algorithms were constructed to solve the governing 
equations. The first algorithm was developed by modifying the unsplit 
finite difference scheme of MacCormack. The second algorithm employed a 
hybrid pseudo-spectral technique in the normal direction to the flow for 
improved resolution of the reacting flow field. The finite difference 
scheme was used in the streamwise direction. It was suggested that more 
attention be given to the development of spectral methods that could be 
more easily applied to high gradient regions like shocks in supersonic 

reacting flows. The case considered with the spectral method was 

shockless and a small degree of damping was applied in the regions of 
high gradients. Several important conclusions were drawn from this 
study, and interested readers should refer to Ref 5. Here, it is 

important to point out that the use of a more complete chemistry model 

rather than the global model in the fluid dynamics equations did not 
result in a set of temporally stiff equations. 



9 


Incorporation of the finite rate chemistry model into the fluid 
dynamics equations can create a set of stiff differential equations. 
The stiffness is due to a disparity in the time scales of the governing 
equations. In the time-accurate solution, after the fast transients 
have decayed and solutions are changing slowly, taking a larger time 
step is necessary for efficiency purposes, but explicit methods still 
require small time steps to maintain stability. An eigenvalue problem 
associated with stiff ordinary differential equations (ODE) has been 
solved to express this point clearly in [29]. The literature related to 
stiff differential equations is not reviewed in this study, but there 
are interesting reviews of this topic available in Refs. 29 and 30. One 
way around the problem is to use a fully implicit method. This method, 
however, requires the inversion of a block multi -diagonal system of 
algebraic equations. The use of a semi-implicit technique, suggested by 
several investigators [28, 31, 32], provides an alternative to the above 
problems. In this technique, the source term which is the cause of the 
stiffness is treated implicitly, and other terms in the governing 
equations are treated explicitly. 

The literature survey indicates that a great deal of effort has 
been directed toward the formulation of radiative heat transfer 
equations. Most applications of these formulations have been restricted 
to nonreacting homogeneous systems. Only a limited number of studies 
are available on applications of radiative energy transfer in combustion 
processes rocket exhaust plume analyses, and fire research. Virtually 
no efforts have been made to investigate the interaction of radiation 
heat transfer in chemically reacting, viscous, subsonic and supersonic 
flows of molecular species. 
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The objective of this study is to investigate the effects of 
radiative heat transfer in chemically reacting supersonic flow in a 
scramjet combustor. The products of hydrogen-air combustion are gases 
such as water vapor and hydroxyl radical. These species are highly 
absorbing and emitting. The presence of such gases makes the study of 
radiative heat transfer in chemically reacting flows an important 
issue. It is essential to employ an accurate and nonosci llatory 
numerical scheme to solve the governing equations involving radiative 
heat transfer. This is because the rate of radiative heat transfer is 
influenced strongly by the temperature and pressure of the medium. 
Consequently, another objective of this study is to explore the 
feasibility of employing the Jameson's four-stage Runge-Kutta time 
stepping scheme [33] for the solution of chemically reacting and 
radiating viscous flows. This scheme has proven to be efficient and 
robust for the solution of the Euler equations. For the steady state 
solution, various techniques can be implemented to accelerate the 
convergence [34-36]. 

A brief discussion of the scramjet engine and the governing 
equations are presented in Chap. 2. Chapter 3 provides the formulation 
for nongray, pseudo-gray, and optically thin radiation models. The grid 
generation technique and solution procedures for the governing equations 
are presented in Chap. 4. Discussion of the results for several 
specific cases are provided in Chap. 5, and specific conclusions and 
recommendations for future studies are provided in Chap. 6. 



Chapter 2 


GENERAL FORMULATION 

A brief discussion is presented on various components of the 
scramjet enoine. Special attention is directed to discussion of the 
basic equations that are applicable in analyzing the flow field in 
different parts of the engine. The relations for the thermodynamic and 
chemistry models are also provided in this section. 

2.1 Physical System and Model 

As mentioned in the introduction, the scramjet engine has been a 
candidate for propelling hypersonic vehicles. In Fig. 2.1, various air- 
breathing and rocket propulsion alternatives are shown [37, 38]. For 
Mach numbers of zero to three, turbojet air-breathing systems have the 
highest performance. Above Mach number of three, turbine inlet 
temperatures constrain performance, and then the ramjet becomes more 
attractive. At about a Mach number of six, the performance of the 

ramjet is greatly reduced. This is due to dissociation of the reaction 
products, which is caused by slowing the supersonic flow to subsonic 
flow through the normal shock that exists in a ramjet. Therefore, it is 
more efficient to allow the engine internal flow to remain at a 
supersonic speed. Thus for Mach numbers of six and beyond, the fixed 
geometry scramjet is clearly superior for propelling a vehicle at 
hypersonic speed. Hydrogen has been selected as the fuel for the 
scramjet due to its capability of cooling the engine and the airframe 
and also because of its high impulse level. 


11 
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The scramjet engine is made up of several identical nodules and it 
is installed underneath the aircraft as shown in Fig. 2.2. Each nodule 
is made up of inlet, combustor and nozzle regions. As part of the 
engine design concept [37], the forebody of the aircraft acts as an 
inlet for precompression and the afterbody as a nozzle for post 
expansion. 

The inlet region starts with the forebody of the vehicle and ends 
up with the minimum cross sectional area of each module. In the first 
part, the air flow is compressed by the oblique shock generated from the 
forebody before it enters the engine. For numerical solution, the flow 
is best represented by the Navier-Stokes equation in the actual inlet 
area of the engine. Using the Euler equation away from the wall and the 
boundary layer equation near the wall region can be complicated by 
oblique shock interaction with the boundary layer. This can cause flow 
separation, which means flow can not be represented accurately by these 
equations. Three-dimensional Navier-Stokes equations have been employed 
by Kumar [1] to investigate the flow field in this region with 
reasonable success. Chitsomboon et al. [4] have employed the 
parabolized Navier-Stokes equations with limited success. 

The combustor region is by far the most complex part of the 
scramjet engine. As a result, a great deal of research is directed 
toward better understanding of the combustor flow field. The flow in 
this region is nearly supersonic, but does have subsonic regions near 
fuel injections (Fig. 2.3). The fluid dynamics become complicated by 
the fuel injection, flame holding, chemistry, radiation and 

turbulence. The flow field in this region is represented by the Navier- 
Stokes equations (including turbulence, chemistry and radiation). In 
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Fig. 2.3 Flow field near an Injector. 
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the farfleld (downstrean of the fuel injection strut, where the flow is 
only supersonic), the flow can be represented by the parabolized Navier- 

Stokes equations [4]. 

The nozzle and subsurface of the afterbody provide about fifty 
percent of the thrust at Mach number six [38]. The flow through the 
nozzle is supersonic. The combustor exit flow consists of multicom- 
ponents of the reacting species, multiple shock, and 3-D viscous 
effects. For engineering accuracy, the flow can be represented by the 
parabolized Navier-Stokes equations. 

2.2 Basic Governing Equations 

The two-dimensional Navier-Stokes and species continuity equations 
are represented in the physical domain by 

i>J + if. + + H = 0 (2.1) 

3t ax ay 

where vectors U, F, G and H are expressed as 
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The viscous stress tensors in the F and G terms of Eq. (2.1) are given 
as 

(2.2a) 


D 1 X 3V\ 9 ou 

T xx ‘ p ^ ax ay^ * u ax 


au 


/3u avv 
T xy “ " wv ay ax' 


(2.2b) 


D . /3U , 9V\ 9 <3V 

T yy = p " X ^x + ay^ u ay 


av 


(2.2c) 


The quantities q cx and q cy in the F and G terms are the components of 

the conduction heat flux and are expressed as 

m 


^cx * - k & - pD J j 1 t<»y»>v 


(2.3a) 


q cy ■ - k £ ' 0° jii [3f J /9y)h J ] 


m 


where 


(2.3b) 


h . - hj ♦ J C dT; T, 


0 K 


It should be noted here that D represents the effective binary diffusion 
coefficent and is used for all species. Assuming that the Lewis number 
(o/D) is unity, Eqs. (2.3) reduce to (see Appendix A) 


= _ m li 
cx Pr ax 

(2.4a) 

- IE 1® 
cy " ‘ Pr ay 

(2.4b) 


where 


e = h - P/p 


The molecular viscosity, is assumed to be temperature dependent and 
it is evaluated from the Sutherland's formula as 

f jJ /2 T °* S 

‘ M 0 M Q ' T + S 


(2.5) 
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where u and T 0 are reference values for individual species and S is the 
Sutherland constant. In this study, reference values were selected for 
pure air because the flow is dominated by nitrogen. The total internal 

energy E in Eq. (2.1) is given by 

2.2 m 

E = P/p + U - ■ o -— + s h. f . (2.6) 

d i = l 

Specific relations are needed for the chemistry and thermodynamic 
models and for the radiative- transport. The chemistry and thermodynamic 
models are discussed briefly in the following sections. The formula- 
tions for radiative transport are presented in Chap. 3. 

Instead of solving the partial differential equations, sometimes it 
is desirable to solve the integral form of these equations using the 
finite volume method. The integral form of the governing equations can 
be expressed as 

JL- fl Udxdy +/ M • n ds + // H dxdy = 0 (2.7) 

9t n an a 

where n is a normal vector pointing outward, n is the region of interest 
and an is the boundary curve. In two-dimensions, the volume has a unit 
depth. The second order tensor M in Eq. (2.7) is defined by 

' M ® F i x + G iy (2.8) 

where F and G are defined in Eq. (2.1). Equation (2.7) can be written 
in Cartesian form as 

|r ff Udxdy + / (Fdy - Gdx) + JJ H dxdy = 0 

dZ n an n 


(2.9) 
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The finite volume formulation is used with the Runge-Kutta time-stepping 
scheme. The details of the formulation and the solution schemes are 

presented in Chap. 4. 

2.3 Chemistry Model 

For the numerical solution of reacting flows, a chemistry model is 
needed to represent the combustion process. The chemistry model used in 
this study is the two-step, global, finite-rate, hydrogen-air combustion 
model developed by Rogers and Chinitz [27]. This chemistry model was 
deduced from a 28 reaction model and is adequate for initial 

temperatures between 1,000 K and 2,000 K and for equivalence ratios 
between 0.2 and 2.0. In the first step, hydrogen and air react and 
produce hydroxyl radical, and in the second step, hydroxyl radical and 
hydrogen react to produce water vapor. The reactions are expressed as 


H 2 + °2 


2 OH 


20H + H 2 



( 2 . 10 ) 


( 2 . 11 ) 


where and represent the forward and reverse reaction rate 
constants respectively. The relations for are obtained from an 
Arhenius equation as 


K f = A.U)T exp(-E./RT) 


( 2 . 12 ) 
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The values of the parameters A..(<|>), N^, and in Eq. (2.12) are 


AjU) = (8.917 4 , + 31.433/4, 

- 28.95)x 10 47 ~ cm 3 /g-mol-S 


Ej = 4865 Cal/mol ; ^ = -10. 

A 2 (<fr) = (2. + 1.333/4, 

- .8334»)x10 64 ~ cm 6 /mol 2 - s 


E 2 = 42500 Cal/mol; N 2 = -13 

where 4 , is the fuel -air ratio. The reverse rate constants can be 
evaluated by 


where 


K 



(2.13) 


Kj = 26.164 exp(-8992/T) 

K 2 = (2.682xl0~ 6 ) (T) exp(69415/T) 


Knowing the reaction rate constant K f and K b the production of species 
can be evaluated from the law of mass action. Consider the general 
chemical reaction 


K 


A li C J 

J 1 K 


f. N 

E C, 


j=l 


J.i J 


i*l ,2 »• • • ,m 


(2.14) 


b i 


where v*. . and v 1 . 1 . represent the stoichiometric coefficents of the 
J > i J » 1 

reactants and products respectively and N is the number of species. The 
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law of mass action states that the rate of change of concentration of 
species j by reaction i is given by [39] 


< c j> = S -\!i) A c j 


ji _ 




N 

n 

j=l 


ji 


(2.15) 


The net rate of production of species j in all reactions is given by 


C = " (Cj), 

J i=l J 


(2.16) 


where m is the number of reactions. Finally, the chemistry source 
terms, on a mass basis, are found by multiplying the molar changes and 
corresponding molecular weight 


w. = t. M . 
J J J 


(2.17) 


By applying the law of mass action to the global model, the chemistry 
source terms of the four species are obtained as 


V 


'OH 


■ - K « c h 2 c o 2 * K b, 4 

(2.18) 

■ 2(K f 2 4 c h, - K b, 4o» 

(2.19) 

2 2 2 


■ E o, • 1/2 Vo 

(2.20) 

2 2 


* - (2 4 * V.> 

(2.21) 


Cj . p fj/Mj 


where 
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2.4 Thermodynamic Model 


The specific heat of individual species, 
linear function of temperature, i.e., 



is assumed to be a 


C_ = a. T + b. (2.22) 

p j J J 

where a^ and b< are constants which are obtained by curve fitting the 
thernochemical data of Ref. 40. The numerical values of these constants 
are given in Table 2.1. The specific heat of the mixture is computed by 
summing specific heats of individual species weighted by species mass 
fraction 


C 

P 


m 


Z C P 
j=l V 3 


f . 
J 


The static enthalpy of the mixture can be expressed as 


(2.23) 


m 

h = I 

j=l 


[h»*j 


T 



dT] fj 


(2.24) 


The total enthalpy can now be evaluated as 

H = h + 0.5 (u 2 + v 2 ) (2.25) 

Combining Eqs. (2.22), (2.24) and (2.25), the total enthalpy is 
expressed as 


m 

H = I 
j=l 


[h j + 


* b j T] 


f j + 


0.5 


(u 2 + 


v 2 ) 


(2.26) 
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where h° is the sensible enthalpy of individual species at a reference 
J 

temperature (T = 0 K). The gas constant for the mixture also is 
evaluated by a mass weighted summation over all the species as 

m 

R = z f. R. (2.27) 

j.i J J 

The equation of state for the mixture of the gases therefore can be 
written as 

P = P RT (2.28) 
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Table 2.1 Numerical Values of Various Constants 


Species 

H°(J/kg) 

a 

b 

°2 

- 271267.025 

0.119845 

947.937 

h 2 o 

-13972530.24 

0.43116 

1857.904 

h 2 

-4200188.095 

2.0596 

12867.46 

OH 

+1772591.157 

0.16564 

1672.813 

No 

-309483.98 

0.10354 

1048.389 


Chapter 3 


RADIATION TRANSPORT MODELS 

In order to include the effects of radiative interaction in a 
physical problem, it is essential to accurately model the absorption- 
emission characteristics of participating species and provide a correct 
formulation of the radiative transfer processes. These are discussed 
briefly in this section. 

3.1 Radiation Absorption Models 

Many models are available in the literature to represent the 
absorption-emission characteristics of molecular species; a review of 
important models is available in [41]. Perhaps the simplest model is 
the gray gas model where the absorption coefficient is assumed to be 
independent of the wavelength. Many nongray models are also available 
in the literature. Both gray and nongray absorption models are 
discussed here briefly. 

3.1.1 Gray Gas Models 

In the gray model, it is assumed that the absorption coefficient is 
independent of the wavelength. This is rarely a physically realistic 
approximation, but it serves as an initial step for studying the effect 
of radiative heat transfer. The absorption coefficient for the gray gas 
is evaluated by employing the Planck mean absorption coefficient defined 
as 


24 
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J k e. (T) d w 


a) 


'bo, 

17(17 


(3.1) 


By assuming that within a band the Planck function does not vary 
significantly with the wave number, and by evaluating its value at the 
band center, the relation for Kp for a single-band gas can be written as 


V T) 

* J da > 


' P a T 4 (y) Au, 


(3.2) 


where u> c represents the band center. For a multiband gaseous system of 
n gases, < p is given by 


n < T > 

tc_ = £ [ 2 J K do) ] 

P k=l a T (y) Au) 


(3.3) 


where represents the band center of the kth band of a particular 

species. For a specific band of a given gas, the integrated band 

intensity is defined as 


1 


S = TT" / < du 

K r . J (1) 

J Au) 


(3.4) 


Substituting Eq. (3.4) into Eq. (3.3), < is expressed as 


= k=i e K 


(T) S k (T) 


(3.5) 


f 3 

/T\ ^ 

e bu>. ^ " exp LC, oj u /U -1 


'2 “V 


where 
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In Eq. (3.5), Pj is species partial pressure, and Cj and C 2 are 
constants [6]. Note that < p is a function of temperature and species 
partial pressure. 

3.1.2 Nongray Gas Models 

Important nongray models available in the literature are as 
follows: 

1. Line Models: 

(a) Lorentz 

(b) Doppler 

(c) Lorentz-Doppler (Voigt) 3 * * * 7 

2. Narrow Band Models: 

(a) Elsasser 

(b) Statistical 

(c) Random-El sasser 

(d) Quasi-Random 

3. Wide Band Models 

(a) Box or Coffin 

(b) Modified box 

(c) Exponential 

(d) Axial 

The relative importance and range of applicability of these models are 
discussed in Ref. 41. In the moderate temperature range (500-5000K), 
use of the wide band models and correlations provide sufficient 
accuracies. These models render significant computational efficiency 
over the line by line or narrow band models. 

The expression for the total band absorptance is given as 

-/ < v U)<is 

A(y> • I [1-e 0 ]dv 

Av 


(3.6a) 
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Differentiation of Eq. (3.6a) gives 

-f r K v (0 

A 1 (r) = / < e 0 dv (3.6b) 

Av v 

< V (C) d 5 

A' ' ( r) = / k e 0 dv (3.6c) 

Av v 

These relations are used to obtain expressions for the total radiative 
flux. 

The radiative flux term usually involves multiple integrals even 
for the simple geometries. As a result, numerical calculations for 
energy transfer becomes very time consuming. Therefore, it is desirable 
to replace the relation for the total band absorptance as given by Eq. 
(3.6a) with a continuous correlation [8, 9, 42]. Several correlations 
are available in the literature for the wide band absorptance. The 

first correlation to satisfy the linear, square-root and logrithnic 
limits of the wide band absorptance was proposed by Edward and Menard 
[43]. The most widely used correlation is the Tien and Lowder 

continuous correlation, and this is given by [42] 

A = in [uf ( 6 ) ( u+ 2 > ~ ^ ) + (3 * 7) 

f(e) = 2.94 [l-exp(-2.6e)] 

0 = B 2 P e 

The form of f(0) was chosen to give agreement with the correlation of 

Edward and Menardo. The Tien and Lowder correlation is employed in this 

study for the nongray gas formulation. 
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3.2 Radiative Flux Equations 

The equations of radiative transport are expressed generally in 
integro-differential form; the integration involves both the frequency 
spectrum and physical coordinates. To overcome the complexity of the 

radiative transport equations, a tangent slab approximation was employed 
in [44, 45]. This approximation treats the gas layer as a one- 
dimensional slab in evaluation of the radiative flux. The multi- 

dimensional equations of radiative heat transfer are formulated for an 
arbitrary geometry, and then an approximate method is used to present 
the formulation for gray and nongray gases. 

3.2.1 Basic Formulation 

The radiative transport equations in the present study are obtained 
only for an absorbing-emitting medium contained within solid walls of 
arbitrary configuration as shown in Fig. 3.1. The general formulation 
of radiative transfer for the gas under the condition of local thermo- 
dynamic equilibrium is given by [12]. 

p* 

-/ w « v (s> is r -I *v ( «> dc 

I V (P) - W e 0 + / W s(P') I bv < p,) e ° dr ' 

0 (3.8) 

The first term on the right hand side represents the contribution from 
the wall to the intensity at P, and the second term, the contribution 
from the intervening gases between P and P w . The origin of the 
coordinate system in Eq. (3.8) is chosen at the point P. 
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The radiative flux at the point P in the direction of t is 


q D = J J I t dv da 

R 71 1 * v 

4tt o 


(3.9) 


Substituting Eq. (3.8) into Eq. (3.9), the radiative flux is expressed 
as 


„ -/ " « V U) d ? 

3 r * f I I V (P W ) <= ° 3 ^ d « 

4-tt 0 


« r w 


-/ * v <5) d « 


*111 %( p ') I bv (f " ) 6 

4tr 0 0 


t dr' dv da (3.10) 


The divergence of the radiative flux is formulated as 


V • = J / t • 7 I dv da 


4 it o 


(3.11a) 


A bean of intensity I (r, a) traveling in the t direction satisfies the 


equation of radiative transfer 


t • 7 I ®y I + < 1 l. ■ & I 

** ' •> " bv v v 


V V V V 


(3.11b) 


where y v is the scattering coefficients and = y v + is referred to 
as the extinction coefficient. For negligible scattering (y v 35 0), Eq. 
(3.11b) reduces to 


t * V I = K (Iu - I ) 
v v bv v 


(3.11c) 
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By combining Eqs. (3.11a) and (3.11c), one obtains 

v • q R = / / < v n bv - V dv dn (3. lid) 

4n o 

Mow substituting Eq. (3.8) into Eq. (3.11b), the divergence of the 
radiative flux is expressed as 

-/ r ” \,< 5 > 

7 • q R ■ 4, /* « v (P) I bv (P) dv - ^ % (P) I v (P w ) e ° 

r' 

r w -/ K v U) 

- / f/\ (P) < v ( pl ) Ibv (p,) 6 0 dr ' dv d " 

4w o o (3.12) 

Equations (3.10) and (3.12) are used to obtain various approximate forms 
for the radiative flux and its divergence. 


3.2.2 


Formulation 


In the previous section, it was observed that the radiative flux 
terms are represented by an integro-differential equation. Solving 
these equations is extremely time consuming because of the complexity of 
integration over space and frequency. Therefore, a pesudo-gray model is 
selected for efficient parametric studies. To express the radiative 
flux for a gray medium, one may assume that k y is independent of the 
frequency. This is rarely a physically realistic approximation; but it 
serves as an initial stepping stone towards nongray analyses. 
Therefore, Eqs. (3.10) and (3.12) for a gray medium are written as 
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and 

v . q R 


-/ W <(K) dc 

4ir <(P) I h (P) - / k(P) I(P w ) e 0 dn 

D 4tt 


r* 

r -/ <U) 

- / / W <(P) k(P') I b (P') e 0 dr' d a 

4 * 0 (3.14) 


To solve the radiative flux terns for the gray medium, one can transforn 
the above equations into the Cartesian coordinates and then apply any of 
the standard integration techniques to evaluate the radiative flux 

terms. It is a lot more convenient and efficient to convert the 

equations to a set of ordinary differential equations (ODE). For the 

present case, differentiating Eq. (3.14) by using the Leibnitz rule 

results in r 

3e -/ W «(?) d 5 

i 2 q p = MR) ♦ / * 2 (P) I(P W ) 6 ° 

R 3r 4u 

r* 

r w -/ <U) d? 

+ / / < W (P) le(P') I b (P') e 0 dr' d a 

4 * o (3.15) 

A substitution of Eq. (3.15) into Eq. (3.13) gives a second order 

nonhomogenious ordinary differential equation as 

- « 2 (P) q R = 4«(P) ^ (3.16) 

dr^ 

It should be pointed out that if the tangent slab approximation is 
employed, then the method of exponential kernal approximation is used to 
convert the equations to a set of ODE's. The result is of the same form 
as Eq. (3.16) but the coefficients are different, i.e.. 
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-t' z < p > q R = 3 ‘ (p) ;fr 

dr 


(3.17) 


Equations (3.16) and (3.17) differ in the coefficients due to the 
enploynent of the exponential kernal approximation in the second 
approach. These equations require two boundary conditions. For 
nonblack diffuse surfaces, the boundary conditions corresponding to Eq. 
(3.17) are given as 


dq c 


11 i l M R i 

(77 ' q R ^ r ) I r=0 ' 1 IT I r=0 = 0 


(3.18a) 


(77 ' 2^ K p q R ^ I r=L + 3 dr I r=L " 0 


(3.18b) 


Detailed derivation of Eqs. (3.17) and (3.18) are given in [44]. 


3.2.3 Nonqray Formulation 

For simplicity in notation, it is assumed that the molecular gas 
has only one band. The following analysis is also applicable to the gas 
having more than one band. Since the Planck function within a wide band 
can usually be approximated by its value at the band center, Eqs. (3.10) 
and (3.12) can be written as 

-/ \,U)dc 

q. - / / I (P w ) [e 0 - 1] J dv dfl + f / I (P w ) t • dv da 

r 4 it Lv V , 4ir Av 

r w -/ P < V U) 

+ f / I. (P 1 ) t dr'da / k (P‘) e 0 dv (3.19) 

i J DvC \ V 

4 7T 0 AV 
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and r 

-I W x v U) d 5 

7 • 3 R = 4, I bve (P)(/ « v dv) - f I vC (P w )[/ <„(P w )d ° <M<*> 

R v Av 4ir Av 

r -/ %U)d5 

+/ / w [I. (P'Jjfcdr'dnJ k (P 1 ) e 0 d v ]dr'da (3.20) 

4tt bvC Av V 

Substituting Eqs. (3.6a-3.6c) into Eqs. (3.19) and (3.20), rearranging, 
one gets 

% ■ l u « P w> 4 d » - > vc ( p w) A(r w > 4 dC! 

r w 

+ J J I b c (P‘) A'(r') t dr'da (3.21) 

4tt 0 ^ 


'•q R 


* 4- I bvC (P)A'(0) 


- J 

4tt 


I vc( P w) A, ( rw)d0 


L 


\vC< P '> A ' 


(r' )dr'dn 

( 3 . 22 ) 


The first tern on the right hand side of Eq. (3.21) represents the 
radiative flux in the absence of participating medium; the second tern 
gives the portion of the wall radiation which is absorbed by the gas, 
and the last tern indicates the emission from the gas which arrives at 
point p. Since I b and A' are independent of n, the first term in Eq. 
(3.22) can be written as 

4. I bvC (P)A'(0) - k I bvc (P)A'(0)da 


■/ I b vc< P > A,(r W > <ln M I bvC (P)CA'(0)-A'(r w )]da 
4tt 4tt 

•I I buC ( p )A'(r )da / / r "l (P)A"(r')dr'd!i (3.23) 

4 it 4ir o 
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A substitution of Eq. (3.23) into Eq. (3.22) results in 

[W P > - I vc< p w»* , ' r «> da 
4tt 

-/ / P " ['bvC< P > * I bvc< p,) ] fl " tr ' )dr ’d« < 3 - 24 ) 

4tt 0 

The first term on the right hand side of Eq. (3.24) represents the net 
exchange of radiation between the gas element at P and the walls; the 
second term represents the net exchange between the gas element of P and 
the other gas elements. 

For a nongray planer system with two parallel walls at different 
temperatures (Fig. 3.2), Eqs. (3.21) and (3.24) reduce to 


« r ■ . x - - J* /* /\ c , 

—oo —oo 0 1 7T r i 


l 1 A -,' 


+ /"/"/ L e <«'• z ' )ft( r 2 ) ^ 

-oo -oo y C 7T l a 


+ \ \ f e (x‘, y', z') £!ip. dx' dy 1 dz' (3.25) 


—00 — OO 0 


nr. 


and 


7-q R ■ CA' (Tj) + A'(r 2 )] e vC (x, y, 2 ) 


• /* " V (x '’ 2,) A ' (r l ) dx ' dz ' 

-00 -00 1 7T r 1 
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- /" /” e vc , (x '’ z '> A ' (r 2> J 5 dx ' dz 

-00 -00 2 7T I r\ 


- \ /“ / L [e c (x. y. 0 - e vC (x '’ y - 2 ')] d ‘ x ‘ dy ' d2 ' 


— 00 —CO 0 


7T r 


(3.26) 


where 


rj = [(x-x') 2 + (y) 2 + (z-z') 2 ] 


1/2 


r ? = [(x-x') 2 + (y-L)‘ + (z-z 1 ) 6 ] 


m2 


1/2 


= [(x-x’) 2 + (y-y') 2 + (z-z 1 ) 2 ] 


1/2 


The y coordinate is nornal to the wall and x and z coordinates are in a 
plane parallel to the wall. For one-dimensional radiative transfer, 
Eqs . (3.25) and (3.26) reduce to 


■ e l - *2 - e vCj A(y) * e vc 2 A(L - y) 


+ / e (y ’ ) A' (y-y 1 )dy' 

o VL 


L 

/ e (y*) A' (y'-y)dy' 
y V (3.27) 


V • \ ■ [« vC (y) - e vCi l A '(y) * [e vC (y> - e^] A'(L-y) 

- \ Iv (y) ' c vc (y,) ! A " < y - y,)dy ' 

0 

- / L lV (y) ■ %C (y) - e vC (y,) ] A " (y '- y)dy 

y 


(3.28) 



Equations (3.27) and (3.28) can be rearranged and written for a multi- 
band gaseous system as (see Appendix B) 
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and 


q R (y) = e, - e 2 * 


N y 


z 

i=l 


{/ F 


IvC • 




(y 1 ) Aj(y-y') dy 1 
(y‘) A’, (y'-y)dy'} 


v*g R (y) 



1-1 

y 


+ / 


0 


f[F lv c (»1 * F 2vC .W1 / 
1 i 

Fj Cy') a ; 1 (y-y') dy' 



♦ J L F 2 (y') Aj' (y'-y) dy' 
0 1 

where 


(3.29) 


(3.30) 


F lvC (y) = e vC (),) ' "vC, 

and 

F 2vC ( *> = e vC (y) - e vC 2 

The equations resulting from employing the tangent slab 

approximation and exponential kernal approximation for evaluating the 

exponetial function are of the same general form as Eqs. (3.29) and 

(3.30). These equations are derived in detail in Ref. 9 and are 

expressed here as 
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q R (y) * ej - e 2 {/* Fj„ c .(y') a; [7 (y-y')l «y' 

- ^ F 2u,c/y’> A i [f ( y'-y ) ] d '} < 3 - 31 > 

and 

7-q R (y) ■§ U F i uCf (y) ♦ ^c/y’l ^ \ d “fl 
♦f 2 t/ F i„c.'y'» A i' [| (y-y')]^' 

4 1-1 o 1 wC i 

* / L F 2 UCj (y,) A i' [| (y'-y'ixy't < 3 - 32 > 

Equations (3.29-3.32) are in proper form for obtaining the nongray 
solutions of molecular species. However, in order to be able to use the 
band model correlations, these equations must be transformed in terms of 
the correlation quantities. For the temperature range considered, the 
radiative process is in the optically thin limit [39]. As a result, 
there is no significant difference between the two approaches. It 
should be noted here that for nongray gases, the divergence of radiative 
flux is used as a source term in the energy equation. This is more 
convenient and also avoids scheme dependency in the computation. The 
correlation quantities and details of transformations are given in Ref. 
20. After the transformation, Eq. (3.32) is written as 
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dq D (u) q n u 3 

■!**„(/ F i„ (u ' )A i' i? (u i • u ; ,] du i 

i =1 i o i 


du 


* I 1 F 2u (u') a:' [§ (U- - u.)] dujl 


U. 1 

1 


+ ? i=i A °i + s (u)] 


(3.33) 


Note that A 1 ' ( u ) expresses the second derivative of A(u) with respect 
to u and 

dq R _ du [ PS ( T )/A ] — 

dy" ~ du dy u/ V du 

By defining ij- = £ and Eq. (3.33) is expressed as 


.2 /I 2 r/ p /.,i \ xii r3 JL 


dq R (y)/dy = 9/4 e A u 2 /l {/ Fj (y 1 ) [? ~T (y " y,) l dy ' 

K i=l u i i o i 


L 3 V 

+ / ( yl ) [l-c 1 (y'-y)] dy'} 


y 2u i 


_ A u _ 

3 " °i °1 




(3.34) 


It is often desirable and convenient to express the proceeding 
equation in terms of A 1 rather than A 1 '. This is accomplished by 
integrating Eq. (3.34) by parts. The detail of the integration by parts 
is given in Ref. 9 and the result is given by 
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dq r „ n V V y deu.(y') 


u 

3 °i 


* -Vt/ -V*; ^ (y - y,,ldy ' 

J i=l o 


L deu.Cy') 3 o. 


-/ -v- 


(3.35) 


Equation (3.35) and the Tien and Lowder correlation given by Eq. (3.7) 
can be used to evaluate the radiative flux. 


3.2.4 Optically Thin Formulation 

In the optically thin limit, the fluid does not absorb any of its 
own emitted radiation; however, it does absorb radiation emitted from 
the boundaries [6-9]. In general, the optical thickness of an 
absorbing-emitting system is defined as 

L 

t = f k dy (3.36) 

Oo) j q a) 

If it is assumed that < is independent of the temperature, then Eq. 

0) 

(3.36) is expressed as 


T 


Ou) 


ic L 
a> 


(3.37) 


A radiating system is considered to be in the optically thin limit when 

T < < 1. 

Ou) 

In the optically thin limit, the expressions for the radiative flux 
and its divergence can be obtained from the general expressions by 
following the procedure outlined in [6]. For the gray gas approximation 
and black-bounding surfaces, Eq. (3.14) is expressed for one-dimensional 
optically thin radiation as 
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dq p 
dy 


2 < [2 e. (y) 

03 D 

03 



(3.38) 


For the pseudo-gray gas nodel , is replaced by <p through use of Eqs. 
(3.1) - (3.5). 

For the case of nongray gas formulation, A(u) = u and A'(u) = 1 in 
the optically thin limit [6, 7]. Consequently, Eq. (3.35) reduces in 
this limit to 



It should be pointed out that it is considerably more efficient to 
solve numerically the system of governing equations for the case of 
optically thin radiation than for the general case. In a preliminary 
study, the optical thickness of the present physical problem were 
calculated. It was found that the optical thickness varied between 
0.003 to 0.4 in the temperature range of 1,000 to 5,000 K and pressure 
between one to three atmospheres. 



Chapter 4 


METHOD OF SOLUTION 

The grid generation technique and solution procedures for the 
governing equations using the unsplit MacCormack [46] technique and 
modified Runge-Kutta time stepping scheme of Jameson [32] are briefly 
discussed in this section. 


4.1 Grid Generation 

Grids are generated using an algebraic grid generation technique 
developed by Smith and Weigel [47]. From the computational point of 
view, it is desirable to have a uniform rectangular grid enclosed in a 
cube, where the exterior of the cube represents the physical boundaries. 
To have such grids, the body-fitted coordinate system is transformed 
linearly from the physical domain (x, y) to the computational domain 
(5, n) as follows: 


X 

II 

X 

'•Pnr 

O 

Lower 

(4.1a) 

y l = y u,o) 

Boundary 

X 2 = x (5,1) 

Upper 

(4.1b) 

y 2 = Y (5,1) 

Boundary 

(5,1) n + X (5,0) 
(5,1) n + Y (5,0) 

(1-n) Between the 

/. x Boundaries 

(1-n) 

(4.1c) 
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where 


0 < 5 < 1 ; 0 < n < 1 


The grid should be concentrated in the regions of high gradients to 
accurately predict the solution. Therefore, more grid points are 
required near the solid boundaries. The concentration of the grid 
in the n-di recti on can be accomplished by 


( 0 y +l) - (B y -1) exp[-C( n -l+o)/(l-a)] ^ ^ 

" = ( 2 a+l ) { 1 +expLlC (n-l+a)/ ( 1 -a) _]} " (4 ’ } 

where 6 +1 

c * tn (/rr) 

y 

If a is equal to zero (a=0), the compression takes place only near the 

lower wall (n= 0 ), and if a is set equal to one half (a=l/ 2 ), the 

compression takes place near both walls. The term e y has a value 
between one and two, and as it gets closer to one, the grid becomes more 
concentrated near the walls. Employing this concentration, Eq. (4.1e) 
is written in terms of n as 

X = X ( 5 , 1 ) n + X ( 5 , 0 ) (1-n) 

y = Y (5,1) n + Y (5,0) (1-n) (4*3) 


where 


0 < n < 1 


It should be noticed that the grid is concentrated in the normal 
direction to capture the boundary layer and kept uniform in the flow 


di recti on. 
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4.2 Solution of the Govening Equations 

4.2.1 Modified MacCormack's Finite-Difference Scheme 

The governing equations, Eqs. (2.1), are expressed in the 
computational domain as 


where 


+ it + + h 

at 35 an 


U = UJ 


F = Fy - Gx 
n n 


G = Gx^ - Fy^ 


= 0 


(4.4) 


H = HJ 


0 = X C y n ’ X n 


Equation (4.4) is discretized temporally and written as 

A 

u 


'■ n+1 = u n - At + iSU. + H n+1 ] 

u '■N _ -* 


(4.5) 


35 3n 

The source term H n+1 must next be linearized. It is expanded in a 
Taylor series in time to give 


H n+1 = H n + At |jr + 0 (At) 2 


(4.6) 


or 


H n+1 = H n + At & ( U n+1 - U n , 


au 


At 


(4.7) 


A substitution of Eq. (4.7) into Eq. (4.5) gives the temporally discrete 
equation in delta form as 
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[I + At -4] AU n+1 = -At [|| + || + H] P (4.8) 

3 U 

A 

where U n+1 - U n is expressed as aU n+1 , i? is the Jacobian matrix of H 

au 

and I is the identity matrix. The components of the Jacobian matrix are 
given in Appendix C. 

Once the temporal discretization used to construct Eq. (4.8) has 
been performed, the resulting system is spatially differenced using the 
unsplit MacCormack predictor-corrector scheme. This results in a 
spatially and temporally discrete, simultaneous system of equations at 
each grid point. Each simultaneous system is solved using the 
Householder technique [48, 49] in combination with the MacCormack 

technique, which is then used to advance the equations in time. The 
modified MacCormack scheme then becomes 


[i * 4 t (4)1j] 


** 4 + £ * fi ]?o < 4 - 9a > 


C ■ °?j * 4 


(4.9b) 


3U 1 J 


(4.10a) 


U 1 ?! 1 = U 1 ?. + 0.5 
ij iJ 


[ A °ij +1 




(4.10b) 


Equations (4.9) and (4.10) are used to advance the solution from 
time n to n+1, and iteration process is continued until a desired 
integration time has been reached. 
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4.2.2 Modified Runqe-Kutta Finite-Volume Technique 

In 1981, Jameson, Schmidt and Turkel [33] proposed an explicit 
finite volume technique using Runge-Kutta time stepping scheme for the 
solution of unsteady Euler equations. The scheme is a modified version 
of the classical four stage Runge-Kutta technique. The method is 
fourth-order accurate in time for linear equations and second-order 
accurate for nonlinear equations. The scheme is second-order accurate 
in space for both linear and nonlinear problems, provided the grid is 
sufficiently smooth. This scheme was extended by Swanson and Turkel 
[35] to solve the thin layer Navier-Stokes equations for transonic flow 
over the airfoil. In the present work, this scheme is being extended to 
solve supersonic chemically reacting and radiating viscous flows. The 
governing equation, Eq. (2.9), in integral form for a region a with 
boundary aa is rewritten here for the convenience as 


i_ // U dxdy + f (F dy - G dx) + // H dxdy = 0 (4.11) 

3t a an a 

The discretization procedure follows the method of lines in decoupling 
the approximation of the spatial and temporal terms. The computational 
domain is divided into quadrilateral cells (Fig. 4.1), and Eq. (4.11) is 
applied to each cell separately to obtain a system of ordinary 
differential equations. The resulting equations are solved by the 
finite volume scheme of Jameson. 

Applying Eq. (4.11) to an arbitrary cell, ABCD, and approximating 
the integrals by the midpoint rule, one obtains 


d 

It 


(A fj U (J ) * L U (j 


+ A. . H. . = 0 
ij 1J 


(4.12) 
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Equation (4.12), with the inclusion of artificial viscosity, is written 
as 


±_ ( A .. 

dt ia ij 


V 


+ L U. . 


• DU .j 


H 1j 


= 0 


where L is the spatial discretization operator, and 


(4.13) 


L U,J = FLX ab * FLX bc * FLX cd * FLX DA (4.14) 

The components of Ujj are now cell averaged quantities; A^- is the cell 
area of ABCD. The vector FLX represents the fluxes through the cell 
sides. For example, FLX AB is written as follows 


FLX AB = F AB aY AB " G AB aX AB 


(4.15) 


where F AB and G AB are viscous and inviscid fluxes through the side AB 
and 


aY AB Y B " Y A ; aX AB X B ‘ X A 
Inviscid fluxes are evaluated from two adjacent cells as 


(F AB^i 


nv 


< F «j * F 


ij-1 


W 2 


(4.16) 


To evaluate the viscous fluxes on the cell sides, one must evaluate u x 
and Uy. These components are evaluated from Green's theorem as follows: 


(u ) AB ■ (u x ), , ,,, ■ J— ! // U dxdy ■ A f u dy (4.17a) 

x AB x a. --i/o n x 1-1/2 an 


i , j-1/2 a 


i,j-l/2 an 


(dy) 


AB 


= ( u y>1. j-1/2 = J 


1 


i, j-1/2 0 


// u v dxdy = - 


1 


’1, j-1/2 30 


J) u dx (4.17b) 
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In discrete form, u y for the cell side AB in Fig. 4.1 is written as 


(“x>AB ' ITT" < U E iV E * U B 4Y B * U F lV F * U A lV A> < 4 ' 17c > 

1 , J “1 /c 


Velocities at the cell corners are evaluated from the four adjacent 
cells, and j_i/2 is evaluated as 


i,j-l/2 


7 <* 


. . + A. . . ) 
i,J i.J-1 


4. 2. 2.1 Boundary Conditions - So far, only the interior cells have 
been considered. The implementations of the boundary cells are consider 
in this section. These include supersonic inflow and outflow, and slip 
and nonslip wall boundary conditions. To determine the inflow 
quantities, the flux vectors will be evaluated based on the free stream 
conditions. For outflow conditions, the quantities are extrapolated 
from the interior points. For the inviscid flow at the solid wall, a 
zero wall flux condition is implied, that is 


^ lS * 0 (4.18) 

n 

where V n is the normal wall velocity and aS is the wall cell surface 
area (Fig. 4.2a). Equation (4.18) is written in Cartesian coordinates 
as 


V n aS = u Ay - v ax 

Thus for inviscid flows, the net fluxes through face 2 (Fig. 4.2a) are 
J(F dy - G dx) = 0 Continuity Eq. 

J(F dy - G dx) = (P A y) X-nom. Eq. 
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J(F dy - G dx) = (-Pax) Y-mon. Eq. 

J(F dy - G dx) = 0 Energy Eq. (4.19) 

Pressure at the wall is evaluated by 

V p i,i -f ay + 0 (iy)2 

The pressure gradient is calculated using the three point differencing 
method. For viscous flow, pressure gradient is set to zero (boundary 
layer assumptions). 

Finally, the no slip boundary conditions for viscous flows must be 

considered. The wall fluxes are divided into inviscid and viscous 

parts; the inviscid parts are evaluated from Eqs. (4.19). For viscous 

flow, the velocity component tangent to the wall must be zero. As the 

result of this constraint, the viscous flux vector is nonzero. For 

example, -^and |i:an be computed on face 1 in Fig. 4.2b as 

ay ay 


9jj U A aX A + U B aX B * U c aX c 
A ABCD 

aT _ T A aX A + t b aX b * T c aX C * T P aX d 
w a abcd 


(4.20a) 


T D = °* 5 ( T i, 


1 + T i+l,l ) 


for an adiabatic wall 


For face 4, the gradient of T and u are evaluated from the ghost cell 


CDEF. 
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and 


3T _ 
3y 

T D = 
T n = 


T c aX c + T F AX F + T e aX e + T d AX 0 

* 7 ^ 


T F = T E = T 1,1 


T F * T w 


(4.20b) 

for an adiabatic wall 

for constant wall temperature 


T E* T i,l +2(T w- T 1,l ) 


Ur AXr ^ Up AX/» 

3u EE C C . ,, = _ ,, 

j » u E u ( 

A i ,1 


ay 


Note that aX is direction dependent (aX c = - aX e ). 

4. 2. 2. 2 Artificial Viscosity - The finite volume scheme (Eq. 
(4.12)), unlike the MacCormack scheme, is not inherently dissipative 
and, therefore, it does allow undamped oscillations with an alternate 
sign at odd and even grid points. To prevent large oscillations caused 
by a discontinuity, some kind of artificial damping must be added to the 
scheme. The original damping proposed by Jameson et al . [33] is a blend 
of second- and fourth-order differencing. The basic idea is to add the 
fourth-order dissipative terms throughout the domain to provide a base 
level of dissipation sufficient to prevent nonlinear instability, but 
not enough to prevent oscillations in the neighorhood of shock waves. 
For the linear problem with central differencing schemes, the 
neighboring points decouple. This odd-even decoupling prevents the 
possibility of driving the residual to machine zero. For the nonlinear 
equation, the values are evaluated at the cell sides before evaluating 
the fluxes. This nonlinearity couples all the neighboring points 
together. However, this coupling is weak and convergence to a steady 
state can be slow [33]. In order to capture the shock waves, additional 
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second-order dissipative terms are added locally by a sensor designated 
to detect the shock waves. In recent years, Jameson has modified the 
original damping to lead to a scheme which will behave locally -like a 
TVD scheme [50, 51]. Both dampings were tested in this study by solving 
the two-dimensional Euler equations. The TVD damping prevented the pre- 
and post-shock oscillations which were observed when the original 
damping was used. The results of the two dampings are compared in the 
next chapter. The original dissipation and the TVD dissipation schemes 
modifications are outlined here. To preserve the conservation form, the 
dissipative terms are generated by dissipative fluxes. The damping term 
in Eq. (4.13) is calculated from 


D i,j (U) = d i+l/2,j d i-l/2,j + d i , j+1/2 ' d i,j-l/2 


(4.21) 


where the dissipative flux &\+\/ 2 , j defined by 

d i+l/2,j = e i+l/2,j R i +1/2 ,j (U i+l,j ' U i,j ) 

‘ e i+l/2,j R i+l/2,j (U i+2,j ' 3 U i+l,j + 3 U i ,j 


U i-l,J^ 

(4.22) 


Here, z\l\ /2 is an £ i+l/2,j are ada P tive coefficients and Rj+i/2,j is a 
coefficient chosen to give the dissipative terms the proper scale. An 

appropriate scale is [50] 


n = 1 r A i+lJ- + 

i+l/2,j 1 At i+i j j At i,j 


(4.23) 


An effective sensor of the presence of a shock wave can be constructed 
by taking the second difference of the pressure. Define 


ij 


P i+l,j " 2P iJ * Pi - 1 »J 
P i + U + ZP i,j + 


(4.24) 
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and set 


v i+l/2,j = max ( v i+2,j * v i+l ,j * v i , j ' v i-l»J 


.) 


The coefficients j an< * e i+l/2 j in ori 9^ na ^ ^ orm were writ ^ en 

as 


£ 


( 2 ) 

i+l/2,j 



i +1/2 




(4.25a) 


and 


where 


e m/2J ■ max[0 -- (k<4) • e fH/2.j )] 


k( 2) = 0.25 and K< 4 ) = 1/256 


(4.25b) 


The nodified coefficients and that have approximately the TVO 

property are written as 

<?• k<2) (4 - 26a) 


and 


e (4) 

*1+1/2 ,J 


max(0., K 


(4) 


- a 


V i+l/2,j ) 


(4.26b) 


where 

k ( 2) = i, k (4) = 1/64, and a = 2 

In a smooth region, is proportional to the square of the mesh width 

and e^ 4 ^ is of order one, therefore, dj + 1 / 2 , j in Eq. (4.22) is of order 
three. In the shock region, is zero, and, therefore, the fourth- 
order damping is cut off to prevent oscillation, and e is of order 

one, so that the scheme behaves locally like a first-order scheme. 
However, this does not effect the global second order accuracy of the 
finite volume scheme [35]. 
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4. 2. 2. 3 Time-Stepping Schene - Equation (4.13) is discretized 
temporally and written as 


U n+1 U n 

A ij ( \t * L U u - 0 U 1j * A ij H u = 0 


(4.27) 


i j At 

Substituting Eq. (4.7) for the source term, Eq. (4.27) is written in 
delta form as 


(i * ^ (f)"] 


AU 


n+1 


= -^[L 

A ij 


u" 


ij 


- D U 




Ai . «■: 


n .i 

ij j 


(4.28) 


To advance Eq. (4.28) in time, the modified four-stage R-K technique in 
combination with the householder technique [46, 47] is employed. At 
time level n, the scheme is written as 

U<°> = u (n) 

[I ♦ it (|g)] iU k * - . k [L l)*' 1 - 00(0) t H?j ] 

^ j 

U K = U° + iU K 


K = 1,4 oij = 1/4, - 1/^ > ~ 1/^ > a 4 ^ 

y(n+l) _ y(4) 

For efficiency purposes, the natural and artificial viscosities are 
evaluated at the first stage and frozen for the remaining stages. 

For time accurate solution, the computational time step, At, must 
satisfy the smallest time scales of the fluid and chemistry, i.e., 
At = min (At f , At h ). If the steady state solution is sought (as in 
this study), it is possible to speed up the convergence by using a 
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larger time scale due to the so-called preconditioning matrix (left hand 
side bracket in Eqs. (4.8) and (4.28)), the purpose of which is to 
normalize the various time scales so that they are of the same order 
[28]. To further speed up the convergence, the solution is advanced in 
time with a time step dictated by the local stability limit. In the 
present work, the local At is based on the Courant-Friedrichs-Lewy (CFL) 
stability limit. Local time stepping allows faster signal propagation, 
and thus faster convergence. 

The radiative flux term is evaluated for both gray and nongray 
gaseous systems. In the nongray gas formulation, the divergence of the 
radiative flux is evaluated using a central differencing scheme and is 
treated as radiative source term in the energy equation. Since the 
radiative flux term is in integro-differential form, unlike the other 
flux terms which are only in a differential form, it is uncoupled and 
treated separately. In the gray gas formulation, Eqs. (3.17) and (3.18) 
are discretized by central differencing, forming a tridiagonal matrix 
(see Appendix D). This tridiagonal matrix can be solved efficiently by 
the Thomas algorithm. The radiative fluxes and chemistry rates are 
evaluated at the first stage and frozen in the remaining stages. 



Chapter 5 


RESULTS AND DISCUSSION 

Based on the theory and the computational procedures described 
previously, two algorithms were developed to solve the two-dimensional 
Navier-Stokes equations for chemically reacting and radiating supersonic 
flows. The performances of two damping schemes are compared by solving 
the Euler equations for supersonic flow through a channel with a ten 
degree compression-expansion ramp (Fig 5.1). Then, the Navier-Stokes 
equations are solved by a finite difference and a finite volume scheme, 
and results are compared with another computational method. Finally, 
the extent of the radiative heat transfer in supersonic chemically 
reacting flows is investigated. 

For simplicity in the rest of this discussion, the original damping 
is referred to as Dampl and the modified TVD version as Damp2. In the 
numerical experiment with Dampl, it was found that the tangential 
component of the second order dissipation (e v ') becomes large at the 
inlet region and this causes the flow to separate near the boundary. 
The excessive damping is because of the pressure jump caused by the 
leading edge shock. This behavior was also observed by Turkel [34] at 
the leading and trailing edges of the airfoil. He suggested multiplying 
the viscosity by (M/M ) 4 . Chen et al . [52] suggested multiplying 
£ ( 2 ) by a linear factor which is zero near the boundary and one in the 
farfield. The coefficients K< 2 ) and K< 4 ) must be readjusted each time 
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by changing the free stream conditions or the grid size. The above 

problems make Dampingl less desirable. On the other hand, Damp2 leads 

( 2 ) 

to a scheme which will behave locally like TVD scheme as long as e 
in Eq. (4.26a) is equal to 1/2 in the neighborhood of a shock wave. 

5.1 Non-Reacting Flows 

The non-reacting flow equations are solved first for inviscid and 
then for viscous flows. The inviscid flows are considered for comparing 
the performance and accuracies of two damping schemes (Dampl and TVD 
Damp) in the finite volume scheme. The results are compared with 
calculated results of the inviscid flow by the finite difference scheme 
of MacCormack. In the remainder of this chapter, wherever finite 
difference and finite volume schemes are referred, the reference is to 
the MacCormack and Jameson schemes, respectively. 

The Euler equations were solved for ideal gas flowing at M^ = 5, 
T = 293 K and P = 1 atm in the channel with a compression- 

expansion ramp (Fig. 5.1). A 51x51 grid with uniform spacing in both 
the flow and normal direction was used to solve the flow. Figures 5.2- 
5.4 show the results for the temperature, pressure, and density 
variations, as a function of x for three locations across the channel 
(lower boundary, center of the channel and upper boundary). The results 
are obtained by employing the finite volume scheme with Dampl. Similar 
results are illustrated in Figs. 5. 5-5. 7 when Damp2 is employed in the 
finite volume scheme. Comparing the results of Figs. 5. 2-5. 4 with the 
results of Figs. 5. 5-5. 7, it is seen that pre- and post-shock oscilla- 
tions are removed by employing Damp2. As mentioned previously , 
employing Damp2, the scheme behaves locally like a TVD scheme as long as 







5.3 Pressure variation with x for inviscid flow by finite volume scheme 









scheme. 
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the coefficient is set equal to 1/2 in the neighborhood of a shock 
wave. The temperature, pressure, and density evaluated by the finite 
difference scheme are plotted at the lower wall, center of the channel, 
and the upper wall in Figs. 5.8-5.10, respectively. The density at the 
lower wall is slightly overpredicted in comparison to ideal gas flowing 
through a ten degree shock; consequently, the temperature is slighly 
underpredi cted. 

For the solution of the viscous flow, two test cases were selected. 
The first case considered in that of a supersonic flow over a flat plate 
and the second case is of a supersonic flow in a channel. The simplest 
way to test the behavior of a Navier-Stokes solver is to solve for the 
flow with a low Reynolds number over a flat plate. In this case, a 
temperature equal to the stagnation temperature was specified at the 
wall. Using the fluid properties specified in Table 5.1 and a uniform 
grid distribution, variations in the temperature and the two velocity 
components were calculated at the exit plane and these are plotted in 
Figs. 5.11-5.13. The profiles show very good agreement with the 
calul ations performed by Carter [53]. The oscillation observed in Figs. 
5.11 and 5.12 are due to the bow shock from the plate leading edge. 

In the second case, the solutions are obtained for the supersonic 
flow in the channel with a compression-expansion ramp. The free stream 
properties considered are = 5, T^ = 293K, and = 0.1 atm. The 
corresponding freestream Reynolds number at the exit plane is about 
1.1x10^. This is a significantly more difficult case compared to the 
first case. To capture the boundary layer, the grids are compressed 
near the boundaries; consequently, the grid aspect ratio becomes very 
large. The non-uniformity of the grid creates significant problems 
for some schemes. For example, the finite volume scheme with Dampl 
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Table 5.1 Flat Plate Test Data 


Properties 

Values 

Dinensi ons 

P 

7.0 

N/m 2 

00 

T 

216 

OK 

00 

u velocity 

882 

m/s 

v velocity 

0.0 

n/s 

^eL 

1000.0 


Pr 

0.72 


C P 

1000.0 

J/kg K 

c v 

714.0 

J/kg K 

L 

0.15 

m 

Gri d 

51x51 
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could accurately predict the pressure profiles through the channel, but 
it is not capable of predicting the temperature or velocity profiles. 
It was found that the Jameson scheme with Dampl is capable of solving 
the viscous flow for a Reynolds number of up to ten thousand. For flows 
with Reynolds numbers on the order of a million and above, the viscous 
properties could not be predicted accurately. Dampl was originally 
developed for the solution of inviscid flows on a grid with an aspect 
ratio of order one. For the solution of viscous flows, it is not 
unusual to have a grid with an aspect ratio on the order of two to 

three. This non-uniformity creates an excessive amount of damping in 
the flow direction and insufficient damping in the normal direction. 
Various investigators [36, 54] have suggested different techniques to 
make the damping uniform in both directions. It was found that by 
substituting At and At instead of At in Eq. (4.23), that the damping 
became more uniform; this could also be achieved by substituting p 

instead of P in Eq. (4.24). Of course, this modification requires the 

readjustment of coefficients, and in Eqs. (4.25). By 

employing Damp2 instead of Dampl, the scheme became nonoscillatory, 

robust, and more accurate. The solution was obtained for Reynolds 
numbers of one million and ten million with two different grids. The 
results showed the same trend as the results of finite difference 
scheme. 

Figures 5.14-5.19 illustrate the pressure and tenperture profiles 
for 51x51 and 101x51 grids and for freestream properties of P^ = .1 atm, 
T = 293K, and M = 5. Figures 5.14 and 5.15 illustrate the pressure 

oo oo 

profiles at the lower wall and the center of the channel, and Fig. 5.16 
represents the pressure at the upper wall. The temperature profiles at 
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Fig. 5.17 Temperature variation with x for viscous flow with two different 
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these locations are plotted in Figs. 5.17-5.19. Reducing the number of 
grid points in the flow direction resulted in higher temperatures and 
lower pressure through the shock within the boundary layer. This is due 
to the excessive amount of artificial viscosity in the streamwise 
direction. For optimum results, grid must be compressed in the 
streamwise direction for the high gradient regions to reduce the grid 
aspect ratio. 

The viscous flow properties were calculated by the finite 
difference and finite volume schemes; the results are compared with the 
calculations performed by Chakravarthy [55] using the Roe scheme. The 
pressure profiles are illustrated in Figs. 5.20-5.22 for the same three 
locations across the channel as in Figs. 5.14-5.19. The calculated 

results predict a pressure jump at the inlet, which is caused by the 
presence of a shock from the leading edge of the channel. Downstream of 
the channel on the lower wall there is an increase in the pressure. A 
close examination of the results shov/ed that due to the interaction of 
leading edge shock with the shock and expansion fans from the 

compression-expansion corners, the velocity vector in that region has an 
angle less than one degree with the boundaries. The interaction of the 
flow with the lower wall creates a series of Mach waves which causes a 
small increase in the pressure. On the upper wall in the shock boundary 
layer interaction region, the pressure predicted by the finite 

difference scheme shows an oscillation. The calculated results compare 
very well with the Roe scheme. Figures 5.23-5.25 illustrate the 
streamwise velocity for the same locations as pressure. At the 
compression corner, the MacCormack scheme shows the flow to be 

separated, while the other two schemes predicted the same trend but the 
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flow was not quite separated. The shock interaction with the upper 
boundary creates a small separation bubble. The flow is fully separated 
at X/L x = 0.7; this is predicted by all three schemes. The streanwise 
velocity illustrated in Figs. 5.24 and 5.25 indicate that the Jameson 
scheme predicts a decrease in velocity before the shock and the 

separation region. This is due to tht excessive amounts of artificial 
viscosity in the neighborhood of shocks when the grid spacing is large 
in streamwise direction. To reduce the excessive amount of streamwise 
dissipation in the shock vicinity, the scaling factor (Ecu (4.23)) must 
be reduced. This can be achieved in two ways, one by reducing the grid 
cell area and the other by increasing the time step in the flow 

direction. Reducing the grid cell area in the shock vicinity requires 
more grid points. However, if there are enough points in the flow 

direction, this can be achieved by compressing the grid in the flow 

direction. For the highly stretched grid, it is recommended that R be 
scaled (Eq. (4.23)), with respect to local At x for the flow direction 
and local At y in the normal direction as opposed to At CFL which was 
used in this study. Also, the results calculated using the finite 
volume scheme predict a faster drop in the velocity near the leading 
edge than the Roe's scheme. 

Temperature profiles are shown in Figs. 5.26-5.28. Temperature is 

the most sensitive and difficult property to predict accurately. In the 

finite difference calculations, temperature is strongly dependent on the 

implication of the boundary conditions. In this work, it is assumed 

that the walls are adiabatic. This could be enforced by either 

3T/ay I = 0 or aH/syl = 0 (Appendix E). The above results were 
!y=0 'y=0 

obtained by enforcing the zero temperature gradient. The same trends 
are predicted by all three schemes. The approximate temperature 
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recovery factor for compressible laminar flow with - 5 and Pr 0.72 
is 5.25 [56]. This value is predicted by the Jameson and Roe schemes 
towards the channel outlet. The above recovery factor can be predicted 
with the finite difference scheme by enforcing the zero enthalpy 
gradient at the boundary. Figures 5.29-5.31 illustrate the temperature 
profiles when zero enthalpy gradient was enforced for the finite 
difference scheme. Temperature increases very rapidly at the leading 
edge due to the flow stagnating from the free stream within one grid 
spacing. 

5.2 Reacting and Radiating Flows 

Before proceeding into the evaluation of the radiative heat 
transfer, the performance of the finite volume and difference schemes is 
compared for chemically reacting flows. Premixed hydrogen-air with an 
equivalence ratio of unity flows through a channel (Fig. 5.32) with a 
ten degree compression corner, and freestream conditions, P^ = 1 atm, 
T = 900°K M = 4. A switch is built into the code to prevent 

chemical reaction for temperatures below 1000°K. The results for O 2 and 
H 2 0 are plotted for two different locations across the channel in Figs. 
5.33 and 5.34. Due to the high temperature in the boundary layer, the 

flow is ignited before the shock, but outside the boundary layer there 

is no chemical reaction in the flow before the shock wave. The finite 
difference scheme shows oscillatory behavior near the shock. The 
results predicted by both schemes are seen to be in a good agreement. 

Radiative flux is a strong function of temperature and pressure 

gradients. The temperature gradient over a grid spacing is greater in 
the normal direction than the streamwise direction. Therefore, it is 
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Fig. 5.30 Centerline temperature variation with x for viscous flow (H, 
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Fig. 5.32 Geometry and flow conditions for chemically 
reacting and radiating flow case. 
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Fig. 5.33 Distribution of HoO and 0 ? mass fraction along the channel 
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reasonable to assume that radiative heat transfer is negligable in flow 
direction. However, in general, this is not the case because the global 
gradient of the temperature has a significant effect on the radiative 
flux. Three different geometries are employed for various parametric 
studies. One is a channel with two parallel plates a distance L apart 

(Fig. 3.2); the other is a channel with a ten degree compression- 

expansion ramp at the lower boundary (Fig. 5.1). The third geometry is 

a channel with a compression corner at the lower boundary (Fig. 5.32). 

For the temperature range expected in the scramjet combustor, the 
important radiating species are OH and H 2 O. The spectral information 
and correlation quantities needed for these species are obtained from 
Ref. 20. Both reacting and nonreacting flows are considered. Results 
for the nonreacting, one-dimensional radiative transfer are presented in 
Figs. 5.35-5.38. Completed information on one-dimensional radiative 
transfer with chemical reaction is provided in Ref. 43. Selected 
results of the one-dimensional radiative transfer analysis are discussed 
here briefly. 


For the parallel plate case (3 cm xlO cm), the inflow conditions 
are P =1 atm, T = 1700K, fl = 3.0, and f H n = 0.5, f 0 = .1 and 

00 00 00 n 2 ^ '“'2 

f„ = 0.4. Results for the radiative flux, as a function of the 
nondimensional location along the flow, are illustrated in Fig. 5.35 for 
various distances from the lower plate. It is noted that the radiation 
flux is approximately zero in the center of the channel (y = 1.5 cm) and 
is significantly higher towards the top and bottom of the plates. This, 
however, would be expected because of the symmetry of the problem and 
the relatively higher temperature near the boundaries. The variations 
in the radiative flux are due to the leading edge shock interaction with 


the boundaries. 




O H 7,05 
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The results for radiative flux are illustrated in Figs. 5.36 and 

5.37 as a function of the nondimensional y-coordinate. For P * 1 atm, 
the results presented in Fig. 5.36 for different water vapor 
concentrations indicate that the radiative interaction increases slowly 
with an increase in the amount of the gas. The results for 50% 1^0 are 
illustrated in Fig. 5.37 for two different pressures (P w = 1 and 3 atm) 
and x-1 ocations (x = 5 and 10 cm). It is noted that the increase in 
pressure has dramatic effects on the radiative interaction. The 
conduction and radiation heat transfer results are compared in Fig. 5.38 
for P = 3 atm and for two different x-1 ocations (x = 5 and 10 cm). The 
results demonstrate that the conduction heat transfer is restricted to 
the region near the boundaries and does not change significantly from 
one x-location to another. The radiative interaction, however, is seen 
to be important everywhere in the channel, and this can have an 
influence on the entire flowfield. The results presented in Figs. 5.36- 

5.38 should be physically symmetric, but due to the predictor-corrector 
procedure used in the MacComack's scheme, they exhibit some unsymmetri- 
cal behavior. 

For the parallel plate geometry, a comparison of the divergence of 
radiatve flux for general (nongray), gray and their optically thin limit 
models is presented in Fig. 5.39 for two different y-locations (y = 0.2 
and 1.5 cm). The inflow conditions are = 1 atm, T^ = 1700K, 
M = 4.3. The gray gas formulation is based on the planck mean 

oo 

absorption coefficient which accounts for the detailed information on 
different molecular bands. As such, this approach is referred to as the 
"pseudo-gray formulation." The magnitude of optical thickness 
calculated for this case (0.003 < t < 0.4) shows that the radiation 
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Fig. 5.39 Divergence of radiative flux along the channel for gray and nongray 
models. 
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regime is in the optically thin limit. The gray, nongray, and optically 
thin limit formulation are employed, and results confirm that the 

radiation regime is in the optically thin limit. For the physical 

conditions of the problem, no significant difference in results is 
observed for the two y-locations. The solution of the gray formulation 
in ODE form proves to be about ten times more efficient than the 

solution of the general formulation on the vector processing computer 
(the gray formulation uses 0.056 CRU's per iteration, while the general 
formulation uses 0.57 CRU’s per iteration). The optically thin 
formulation is slightly more efficient than the gray formulation. All 
other results presented in this study have been obtained by using the 

pseudo gray gas formulation. 

To investigate the effects of radiative energy tranfer in 
chemically reacting flows, premixed hydrogen and air with an equivalence 
ratio of unity was selected. Specific results are obtained for one- and 
two-dimensional radiative transfer with the physical geometry of Figs. 
5.1 and 5.32, respectively. For the one-dimensional radiative transfer, 
the inlet conditions considered are = 1 atm, T^ = 1700, = 4.5. 

For the physical conditions of the problem, the radiation participating 
species produced due to the chemical reaction essentially are OH and 
H 2 O. The radiative interaction is started at about X/L x = 0.20 to make 
sure that there are significant amounts of OH and H 2 O produced by the 
reaction for active participation. This restriction was removed later 
and radiative interactions start at the inlet for the remaining 
results. The pressure contours for the flow without and with chemical 
reaction are shown in Figs. 5.40 and 5.41. A comparison of the pressure 
contours shows that the shock angle has increased in the case of the 
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reacting flow. This is due to a relatively thicker boundary layer and 
changes in the thermophysical properties of the mixture. 

The variations in species concentration along the channel are 
illustrated in Fig. 5.42 for chemically reacting, and chemically and 
radiating flows. It is found that due to the radiative interaction, the 
concentration of OH increases by about five percent and the concentra- 
tion of H 2 0 decreases by the same amount. It should be noted that the 
radiative interaction has no significant effect on 0 2 and H 2 . The 
effect of radiatve heat transfer on the temperature is almost 

negl igible. 

The effects of two-dimensional radiative transfer in chemically 
reacting flows was investigated for the physical problem of Fig. 5.32. 
Premixed hydrogen and air with an equivalence ratio of unity, and inlet 
conditions of = 1 atm, T^ = 900 K, = 4.0 were selected. 

Figures 5.43-5.46 are the contour plots of the reactants and products. 
As expected, the destruction and production of the species occurred only 
in the boundary layers and after the shock, where the temperature is 
greater than 1000 °K. Figures 5.47 and 5.48 are contour plots of the 
pressure without and with chemical reaction. It should be observed that 
the compression shock is curved when chemical reaction takes place. A 
series of shock waves is created due to the pressure increase caused by 
a sudden heat release from chemical reaction. Interaction of these 
shock waves with compression shocks increases the shock strength, and 
this causes the shock to curve. 

Figure 5.49 illustrates the variation of nondimensional y-radiative 
flux for several locations across the channel. The radiative flux in 
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Fig. 5.42 Distribution of species mass fraction along the channel for 
reacting, and reacting and radiating flows. 
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Fig. 5.43 O 2 contours for reacting and radiating flow in a channel with 
compression corner. 
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Pressure contours for nonreacting and nonradiating flow in 
channel with a compression corner. 
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Fig. 5.48 Pressure contours for reacting and radiating flow in a channel with 
a compression corner. 
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the y-direction is neglegible before the shock; this is due to the 
almost symmetrical condition and lower pressure and temperature. The 
radiative flux substantially increases after the shock due to the 
increase in participating species, pressure, and temperature. Radiative 
flux decreases away from the lower boundary due to the reduction in 
temperature. The radiative flux in a streamwise direction is 
illustrated in Fig. 5.50. It peaks at a short distance from the 

boundary and then gradually decreases toward the center of the 
channel. The peak is due to the pressure gradient caused by the leading 
edge bow-shock. The radiative flux (q Rx ) remains constant in the region 
with no chemical reaction, and gradually decreases in the region with 
chemical reaction. The reduction in q Rx is caused by cancellation of 
fluxes in positive and negative directions. 

Figures 5.51 and 5.52 illustrate the radiative flux for the same 
inlet conditions as in Fig. 5.32 but for the 15° compression corner. 
Comparing the results of Figs. 5.49 and 5.51, it is noted that, due to 
the increase of shock strength, the radiative flux in y-direction has a 
steep gradient through the shock. After the shock, it gradually 
decreases as some of the radiative flux from lower boundary gets 
canceled from the radiative flux of the upper boundary layer. Figure 
5.52 is the illustration of the streamwise flux along the channel for 
several locations across the channel. Comparing the results with the 
results of Fig. 5.50, it is noted that the radiative flux increases 
considerably with increasing the shock angle. 

Figures. 5.53 and 5.54 illustrate similar results as presented in 
Figs. 5.49 and 5.50 for the same freestream conditions but for = 6. 

the Mach number the shock angle decreases. As a result, 
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Fig. 5.50 Variation of streamwise radiative flux with 
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compression corner (M = 4). 
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Variation of streamwise radiative flux with x in a channe 
compression corner (M = 6). 
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there is no chemical reaction at y = 0.5 cm. Therefore, streamwise 
radiative flux at this location reduces to zero. It is noted that the 
radiative flux increases substantially by increasing the Mach number. 

Figure 5.55 illustrates the variations in species mass fraction for 
two locations (y = 0.02 and 0.2 cm) across the channel. At y = 0.02 cm, 
where the value of total energy is reduced by the q Rx (q Rx direction is 
from outlet to inlet), less H^O is produced. After the shock, more H 2 0 
has been produced, for both locations because the total energy is 
increased by q Ry . It should also be noted that at locations where the 
total energy is reduced less of the reactant (0 2 ) is used (Fig. 5.56). 
The temperature variation for chemically reacting, and reacting and 
radiating flows along the channel are illustrated for two locations (y = 
0.02 and 0.2 cm) in Fig. 5.57. At location y = 0.02, where the value of 
q Rx is greatest in comparison to the other locations, the temperature 
has been reduced. At y = 0.2, there is no change in the temperature 
before the shock because there are no paticipating species. After the 
shock, temperature has slightly increased over the nonradiating case; 
this is due to the contribution from q Ry . 




Fig. 5.55 Distribution of HoO mass fraction along the channel for reacting 
and radiating, and reacting flow. 
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Fig. 5.57 Temperature variation with x for reacting and radiating, and 
reacting flow. 
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CONCLUSIONS 

Based on the theory and computational procedure that have been 
described, two algorithms were developed to solve the continuity, 
momentum, energy, and species continuity equations for nonreacting, 
reacting, radiating and reacting supersonic channel flows. The finite 
volume algorithm with the appropriate damping scheme proved to be 
accurate, non-osci 11 atory and robust for reacting and nonreacting, 
viscous and invsicid flows. The supersonic viscous and inviscid flows 
in a channel with ten degree compression-expansion ramps at the lower 
wall were solved by both the finite-volume and finite-difference 
schemes. The results of viscous flow were compared with the results of 
the upwind scheme of Roe. The finite-volume scheme required more grid 
points in the boundary layer and shocks than did the finite difference 
scheme. There is an excessive amount of artificial viscosity due to the 
highly stretched grid in the boundary layer, near high gradient regions 
(like compression or expansion corners). However, this problem can be 
alleviated by increasing the grid points in the high gradient regions in 
the flow directions. 

The supersonic pre-mixed hydrogen-air flow with equivalence ratio 
of unity was solved by both finite difference and finite volume schemes 
in a channel with a ten degree compression corner. The flow was ignited 
by the shock wave from the compression corner. The results of finite 
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volume scheme compared very well with the results of the finite 
difference scheme. Finite volume scheme did not show the oscillation 
which was observed in the finite difference solution. 

In the hypersonic propulsion system, the temperature ranges from 
one to five thousand degrees Kelvin. In this range, various nonsym- 
metric molecules become highly radiative participating. One-dimensional 
radiative flux was included in the energy equation for the solution of 
nonreacting supersonic flows between parallel plates. It was concluded 
that the radiative flux increases with the increase of pressure, 
temperature, and participating species. In the case of flow without 
chemical reaction, most of the energy transferred was by convection in 
the direction of flow. As a result, the radiative interaction did not 
affect the flow field significantly. 

Finally, two-dimensional radiative interaction was investigated for 
supersonic chemically reacting flow in a channel with a compression 
corner. Some important results were obtained by considering the 
radiative flux in both directions. The results revealed that radiation 
can have a significant influence on the entire flow field; however, the 
influence is stronger in the boundary layers. It was found that due to 
the temperature increase by chemical reaction, radiatve transfer in the 
streamwise direction can have a significant effect on the flow field. 
Radiative heat transfer is a strong function of the pressure 
pathlength. By increasing the dimensions from the model geometry to 
physical geometry, pathlength and, therefore, radiative heat transfer 
will increase. It is concluded that radiative heat transfer in the 
hypersonic propulsion can have a significant effect on the entire flow 
field. It was also found that the numerical scheme based on the pseudo- 
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gray gas formulation for the radiative flux is highly efficient as 
compared to the scheme based on the nongray formulation, especially for 
the vector processing computers (over ten times). It is concluded that 
for the supersonic chemically reacting flows the radiation is in 
optically thin limit for the physical conditions considered in this 
study. 

It is suggested that future study include grid adaption to the 
finite volume scheme for better shock resolution. For steady state 
solution, it is suggested that the multi grid and residual smoothing be 
used to speed up the convergence. The effect of radiative heat transfer 
in non-premixed chemically reacting flow should also be investigated. 
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APPENOICES 



APPENDIX A 


DERIVATION OF CONDUCTION HEAT FLUX TERMS 

To simplify Eqs. (2.3) to Eqs. (2.4), the Lewis number is assumed 
to be unity. This simplification is carried out in detail for Eq. 
(2.3b) and the same is applied to Eq. (2.3a). Using the expressions for 
the thermal diffusivity (a) and Lewis number (L e ), Eq. (2.3b) can be 
expressed as 


a = 


Wp 


«/D 


n 


q =- p D(L C - + I 
H cy M e p ay 


i-1 


3 f 

— h.) 
3y i 


(A.l) 


Defining the binary diffusion coefficient D in terms of the Prandtl and 
Lewis number Eq. (A.l) can be expressed as 


Pr = v/a 


_ a _ v/Pr _ 

" I U 7FFT 

e e v ■ 
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cy 



m 
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i=l 


3f.. 

— M 
3.y i j 


(A. 2) 


where 
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The static enthalpy of the mixture is given by the relation 


n 

h = E 
1*1 



(A. 3) 


It should be noted that n is a dummy variable employed to evaluate the 
sensible enthalpy. Using the Leibnitz formula, Eq. (A. 3) is 
differentiated to obtain 



3C n (n) 


U f 


ay 


+ f i c Pi (T) $ 


(A. 4) 


The coefficients of the first differential on the right hand side is 
equal to 1^ and the second and third terms are identical to zero, 
therefore, Eq. (A. 4) reduces to 


ah. _ 
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af< 
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ll 

ay 


] 


or 


ill - 
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(A. 5) 


Substituting Eq. (A. 5) into Eq. (A. 2), q cy is expressed as 

_ u ah 
q cy ~ ~ Tr ay 


or 


q - -¥■- 

M cy Pr ay 
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APPENDIX B 

RADIATIVE FLUX EQUATIONS FOR NONGRAY GASES 

In this appendix, Eqs. (3.27) and (3.28) are rearranged and written 
as they are expressed in [12]. The first integral in Eq. (3.27) can be 
written as 


/ e„ c (y')A'(y-y')dy’ » / ? [e (y') - e ] A'(y-y')dy' 

n 0 1 


+ e /' A'(y-y') dy' 


'vc 


1 0 


= f [e vC (y') - v,] fl,(y - y,,dy ' ■ e vc/ (y - y,) lo 

o 1 1 


- J [« vC (y,) - e vC J fl ' (y - y,)dy ' • e vc A(0) 
0 1 1 


* e vCj A(y) 


(B.l) 


The second integral in Eq. (3.27) can be written as 


/ L e vC (y')A'(y’-y)dy' = J Cc vC (y 1 ) - A'(y'-y)dy' 

y y 


+ / e A'(y'-y)dy' 
y vC 2 
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Upon substituting Eqs. 
be written as 


t C«„c (r’)-e oy A(/-y)|J 

y 2 £ 

/ Ce vc (y') -e^'ly'-y'Jdy' - « vC ,Aa-y) 

- e A (B.2) 

v c ^ 

(B.l) and (B.2) into Eq. (3.27), the result can 


\ = e l • e 2 ' 6 vC My) + e vC, A(L ' y) + / [e ..<- (y ’ ) ‘ e uC . ]A ' (y - y,) 


'1 


0 vC vC l 


+ e „r fi (y) - / Ce _(y') - e ] A'(y l -y)dy t - e A(L-y) (B.3) 


-vci- w ' J y “'vc w ' 'vc 2 


Equation (B.3) can be sinplified to 


v ■ e i - e 2 * / fv (y,) - \c 3 A '(y-y') d y' - / Ce vc (y') 
1 o i y 


- e vC ]A'(y'-y) dy 


(B.4) 


Equation (B.4) is the same as Eq. (3.29) if it is written for multi-band 
gaseous system. The first integral in Eq. (3.28) can be written as 

y y 

f [ 

0 


J [e (y) - e vC (y , )]A' l (y-y')dy' = / [e - e vC (y')] A"(y-y')dy' 

O 0 1 


-I [e.„ - e.. r (y)]A ,, (y-y')dy' 


vC 

0 1 
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■ A« vC ‘ e vC (y'W'(y-y')dy' * C% c , * e vc (y):l A ' (y - y ' ) lo 

- / [e vc - e (y 1 )] A"(y-y')dy' * C* vc - V(y)]A'(0) 

0 1 * 

- [e v Cl " e vc (y)] A ’ (y) (B,5) 

The second integral in Eq. (3.29) can be written as 

l Ce vC <y> - e vC (y' )]A' ' (y'-y)dy' = / [e^ - e^ty'JA' '(y’-yldy’ 
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-/ Ce if - e _(y)]A'(y-y’)dy 

y VC 2 VC 
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Substituting Eqs. (B.5) and (B.6) into Eq. (3.28) can be written as 


vq r - Ce uc (y) - e vCj ]A'(y) * Ce vC (y) - V^'O-y) 


- / > vCl - e^y'Wty-yW 

0 1 


- e„ r (y)]A'(0) + [e vC - e vC (y)]A' (y) 


'vc^ vC' 


/ C*, - • ,(y , )]A"(y , -y)dy' 

y VC 2 V 


+ tV 2 - e vC^^(^>-Kc 2 - e vC^>'( 0 ) 


Equation (B.7) can be simplified and written as 


(B.7) 


vq r » C(e vc (y) - e vC ^) + ( e vC (y) • e vc 2 )A ' (0) 


+ / ^ e vc (y,) " e v c, ]A " (y - y,)dy ’ + /. [e vC (y,) - e vC , ]A,,(y, - y)dy ’ 


(B.8) 


Equation (B.8) is exactly like Eq. (3.30) if it is written for a multi- 
band system. 
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APPENDIX C 


COMPONENTS OF JACOBIAN MATRIX 

In this appendix, the species matrix in the left hand side bracket 
of Eqs . (4.8) and (4.14) is evaluated. The source term is function of 

gH 

density, temperature and various species. In evaluating jjj, the density 
and temperature dependency is neglected for computational efficiency. 
Eklund et al . [57] evaluated the species Jacobian matrix with and 

without density and temperature dependency and found no difference in 
the steady state solutions. . The components of Jacobian matrix are as 
fol 1 ows: 
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APPENDIX D 

RADIATIVE FLUX EVALUATION FOR PSEUDO GRAY MODEL 

This appendix shows the discretization and evaluation of the 
radiative flux for the y-direction. The sane thing can be applied in 
the streanwise direction. Equations (3.17) and (3.18) in the y- 
direction are written as follows: 


d ^ - | K 2 (y) , r (y> - 3K(y) ^ 


dy 


! dq p (y) 
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K(y) - j)<l r (y) 
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(D.l) 

(D.2) 

(D.3) 


The above equations are discretized by central differencing. The second 
derivative of q r in the physical domain is discretized as 
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j ' 1 6 j 
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Equations (D.1-D.3) in discrete forms are written as 
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[ K l( ei " 2 ) + 3AyJ q l " 3Ay 1 q 2 
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The above can be written in matrix form as 
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The tri-diagonal matrix on the left hand side of Eq. (D.8) is solved 
efficiently by the Thomas algorithm. 
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APPENDIX E 


BOUNDARY CONDITIONS FOR ADIABATIC WALL TEMPERATURES 


The adiabatic wall condition is given as 


3T 

ay 




= o 


Substituting for the total enthalpy, its gradient is written as 


(E.l) 


3H 

ay 



|y (C p T + (u 2 * v 2 )/2) 



(E.2) 


The velocity terns (u, v) are equal to zero at the boundary for the 
viscous fluid. Assuming specific heat is a constant Eq. (E.2) reduce to 


3H 

ay 


y=0 


3T 

ay 


y=0 


(E.3) 


Specific heat is a constant, therefore, temperature gradient must be 
equal to zero. For Pr = 1 total enthalpy remain constant, while for 
variable Prandtl number total enthalpy changes, the gradient remains 
very small. Therefore, at the solid boundary, one can enforce the zero 
enthalpy gradient by setting 

H w ■ H, (E.4) 

Substituting the relation for total enthalpy Eq. (E.4) is written as 

C p T w ■ C p T 1 + < u ( + v l )/2 (E-5) 
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Therefore, the adiabatic wall temperature can be approximated by 


* h * < U 1 * C p ' 



